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Abstract 


The  parameters  of  two  pairs  of  potentials  that  describe  argon  over  its  entire  liquid  phase  at  a 
fixed  pressure  were  optimized  through  a  novel  application  of  constant  temperature  and  pressure 
molecular  dynamics  (NPT-MD)  and  Monte  Carlo  (NPT- MC)  computer  simulations.  The  forms 
of  these  potentials  were  those  of  a  modified  Lennard-Jones  potential  and  a  Lennard-Jones 
potential  (Lennard-Jones,  J.  E.  Physica,  Vol.  4.  p.  941,  1937).  The  optimized  potential 
determined  using  NPT-MD  simulations  reproduced  experimental  densities,  internal  energies,  and 
enthalpies  with  an  error  less  than  1  %  over  most  of  the  liquid  range  and  yielded  self-diffusion 
coefficients  that  were  in  excellent  agreement  with  the  experiment.  The  results  using  the  potential 
determined  by  NPT- MC  simulations  were  in  almost  as  good  agreement  with  deviations  from 
experiments  of  no  more  than  5.89%  for  temperatures  up  to  vaporization.  Additionally,  molar 
volumes  predicted  using  this  potential  at  pressures  in  the  100-600  atm  range  and  over 
temperatures  in  the  100-140  K  range  were  within  0.83%  of  experimental  values.  These  results 
show  that  when  properly  parameterized,  Lennard-Jones-like  potentials  could  describe  a  system 
well  over  a  large  temperature  range.  Further,  the  method  introduced  was  easy  to  implement  and 
was  independent  of  the  form  of  the  interaction  potential  used. 
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1.  Introduction 


One  of  the  major  goals  in  any  computer  simulation  is  to  accurately  predict  the  properties  of 
the  system.  This  can  only  be  accomplished  with  the  proper  choice  of  interaction  potential.  The 
interaction  potential  is  a  fundamental  feature  of  a  system.  All  thermodynamic  and  transport 
quantities  will  be  determined  by  its  form.  Unfortunately,  analytical  functions,  such  as  the 
hard-sphere  potential,  square-well  potential  [1]  and  the  Lennard- Jones  potential  [2],  provide  only 
an  approximate  description  of  the  underlying  interaction  potential.  The  Lennard-Jones  potential 
has  played  an  important  role  in  the  theory  of  classical  fluids.  Its  simple  form  makes  it  desirable 
for  computer  simulations  and  theoretical  calculations  alike.  It  has  been  extensively  used  in 
computer  simulations  of  liquids  [3-5],  glasses  [6-12],  and  phase  coexistence  [13-26],  to  name 
only  a  few.  It  has  also  been  widely  used  as  a  reference  fluid  in  perturbation  treatments  for  more 
complex  fluids  [27],  Further,  several  attempts  have  been  made  to  obtain  an  equation  of  state;  the 
most  notable  being  by  Nicolas  et  al.  [28]  and  another,  and  more  recently  by  Johnson  et  al.  [29]. 

The  Lennard-Jones  potential  requires  only  two  parameters  in  its  description:  £,  which  is  the 
depth  of  the  potential  well,  and  cr,  which  is  the  effective  particle  diameter.  Often,  there  are  many 
parameter  pairs  that  will  predict  with  equal  accuracy  a  thermodynamic  property  at  a  given  phase 
point.  However,  the  potentials  using  these  different  parameter  sets  might  produce  widely 
varying  predictions  at  other  phase  points.  In  principle,  there  should  be  a  single  parameter  set  that 
will  best  describe  the  system  at  every  phase  point  within  the  limitations  of  the  model.  In  this 
study,  this  set  is  referred  to  as  the  “common  set.”  This  report  will  describe  the  method  of  finding 
the  common  set  for  a  modified  Lennard-Jones  potential  and  a  Lennard-Jones  potential  using 
either  molecular  dynamics  or  Monte  Carlo  simulations.  It  will  be  shown  that  a  potential  using 
the  common  set  reproduces  features  of  the  system  over  the  entire  range  studied.  The  Lennard- 
Jones  form  of  the  interaction  potential  can  be  regarded  as  one  member  of  any  number  of 
interaction  potentials  that  could  describe  the  system  properly  over  the  entire  liquid  phase  and, 
therefore,  is  not  unique  for  parameterization  using  this  method.  In  principle,  this  approach  can 
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be  used  to  optimize  any  functional  form  of  the  interaction  potential.  The  method  will  be 
illustrated  using  liquid  argon. 


2.  Simulation  Method 

The  method  presented  can  be  tested  by  performing  either  Monte  Carlo  (MC)  or  molecular 
dynamics  (MD)  simulations  under  constant  temperature  and  pressure. 

2.1  Molecular  Dynamics  (MD).  In  MD  simulations  [30,  31],  driving  a  system  into  a 
nonequilibrium  steady  state  by  coupling  it  to  the  appropriate  external  fields  is  well  described 
[32],  These  nonequilibrium  molecular  dynamics  (NEMD)  simulations  allow  for  faster  and,  in 
some  cases,  more  accurate  determination  of  certain  transport  quantities  than  equilibrium 
simulations.  A  currently  unexplored  concept  is  presented  here,  where  a  nonequilibrated  liquid 
system  at  constant  temperature,  pressure,  and  particle  number  is  driven  to  a  steady  state  under 
the  constraint  of  achieving  a  particular  bulk  property.  For  this  study,  the  bulk  density  was 
chosen  to  be  the  constraining  “external  field.”  Here,  the  external  field  will  continuously  “push” 
the  system  towards  the  desired  bulk  density.  For  these  purposes,  the  desired  density  is  the 
experimental  density  of  liquid  argon. 

Consider  a  nonequilibrium  system  of  N  particles  at  constant  temperature,  T,  and  pressure,  P, 
interacting  through  a  modified  Lennard-Jones  potential, 


v(r)  =  4s 


f  /  \6 

£|  'a' 

^  r  J 


\rj 


+  Ar2  +  B, 


(1) 


where  A  and  B  are  constant  and  chosen  such  that  both  the  potential  and  the  force  are  continuous 
at  the  potential  cutoff,  rcut.  Modifying  the  potential  in  this  way  ensures  that  the  energy  is 
conserved. 
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The  bulk  density  constraint  is  imposed  by  starting  the  system  with  initial  values  of  s  and  o 
and  allowing  s  to  vary  in  time,  t,  according  to 

€(t)  =  €0+X]TAp(t'),  (2) 

t'=i 

where  X  is  the  coupling  strength,  so  is  the  initial  value  (seed)  for  s  and  the  difference  between  the 
desired  and  instantaneous  bulk  densities  is  Ap(t’)  e  pdesired  -  p(t ).  The  effect  is  that  s  will  be 
made  larger,  thereby  increasing  the  attraction  between  particles,  if  the  system  is  “too  dilute.” 
Similarly,  8  will  be  made  smaller  (thus  decreasing  the  attraction)  if  the  system  is  “too  dense.”  A 
value  for  X  that  is  too  large  will  cause  the  system  to  dramatically  overshoot  the  desired  density, 
thus  requiring  a  longer  amount  of  time  for  the  system  to  relax  under  the  “external  field.”  The 
system  momentum  will  be  conserved,  as  in  a  usual  equilibrium  simulation,  since  the  external 
field  is  actually  coupled  into  the  pair  potential  itself.  Failure  to  conserve  momentum  may  be 
further  indication  that  the  value  of  X  is  too  large.  If  X  is  too  small,  the  external  field  will  be 
unable  to  compete  with  the  fluctuations  inherent  in  the  system  and,  therefore,  the  system  will 
never  approach  the  desired  bulk  density. 

Eventually,  the  instantaneous  densities  of  the  system  will  oscillate  around  the  desired  bulk 
density.  It  is  necessary  that  the  system  oscillate  several  times  around  the  desired  bulk  density 
before  8  is  chosen.  This  ensures  that  the  s  that  is  chosen  will  yield  the  desired  bulk  density,  hi 
the  case  of  liquid  argon,  it  was  found  that  the  integration  of  25,000  time  steps  (1  time  step 
=  5.1  fs)  was  needed  for  the  system  to  approach  an  s  that  gave  the  desired  density.  After  this 
time,  an  8  that  reproduced  the  desired  bulk  density  (for  a  given  time  step)  to  within  a  certain 
amount  was  chosen.  This  amount  can  be  defined  with  any  precision  and  was  chosen  here  to  be 
5x10^  particles  A-3.  The  ability  of  the  system  to  come  within  this  established  amount  depends 
on  the  natural  fluctuations  of  the  system. 
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Repeating  this  procedure  over  a  range  of  c’s  for  a  given  temperature  and  pressure  will  result 
in  a  set  of  several  combinations  of  s  and  a’  that  will  reproduce  the  desired  density  for  this  phase 
point.  The  common  set,  which  refers  to  the  parameter  set  that  reproduces  the  density  over  all 
temperatures  in  the  liquid  phase  at  a  given  pressure,  requires  generating  ensemble  averages  of  s’s 
for  the  various  a’s  at  several  temperatures  over  the  range  to  be  studied.  With  the  exception  of 
simulations  for  which  rcut  =  6.00a,  each  ensemble  consisted  of  ten  trajectories;  the  remaining 
ensembles  consisted  of  five  trajectories.  Once  these  parameter  sets  have  been  generated,  the 
average  value  of  s,  <  s  > ,  for  each  a  is  calculated  and  curves  of  <  s  >  vs.  a  for  each  temperature 
are  drawn.  The  a  of  the  common  set  corresponds  to  the  point  where  all  the  curves  intersect 
(within  statistical  error).  The  s  of  the  common  set  is  the  average  value  of  <  £  >  determined  for 
all  the  temperatures. 

The  common  set  obtained  for  the  modified  Lennard-Jones  function  in  equation  1  is 
dependent  on  the  value  of  rcut.  As  rcut  becomes  large,  however,  the  changes  in  the  common  set 
become  less  pronounced,  as  will  be  shown.  In  order  to  examine  this  dependence  on  the  value  of 
rcut  and  satisfy  the  minimum  image  convention  imposed  in  this  simulation  it  was  necessary  to 
increase  the  system  size  as  rcut  increased.  The  systems  consisted  of  500;  864;  1,024;  and 
2,000  particles  for  rcut’s  of  2.50-3.50a  (0.025a  intervals)  4.00a,  5.00a,  and  6.00a,  respectively. 
Particle  masses  were  set  to  39.95  amu,  and  all  simulations  were  run  at  a  constant  pressure  of 
40  bar.  The  pressure  was  held  constant  using  Andersen’s  method  [33],  making  it  necessary  to 
choose  a  piston  mass.  A  larger  mass  will  result  in  better  momentum  conservation  than  a  smaller 
one,  but  the  system  will  explore  volume  space  more  slowly.  In  an  effort  to  satisfactorily 
conserve  momentum  and  sample  volume  space  in  a  reasonable  amount  of  time,  the  mass  was 
chosen  to  be  2  x  10-3  amu  A for  the  systems  consisting  of  500  and  864  particles,  and  5  x  10"4 
amu  A"4  for  the  larger  systems.  Finally,  the  temperature  was  held  constant  through  a  simple 
scaling  of  velocities,  and  the  equations  of  motion  were  integrated  using  a  modified 
velocity- Verlet  algorithm  [30,  34,  35],  which  eliminated  the  need  for  scaled  coordinates.  The 
coupling  strength,  X,  was  chosen  to  be  1  x  10”4  eV  A3  particle-1. 
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Average  values  for  s  were  determined  using  the  procedure  previously  described  and 
equation  1  with  the  different  values  of  rcut.  Each  ensemble,  except  for  the  series  with 
rcut  =  6.00a,  consisted  of  10  trajectories;  the  ensemble  for  rcut  =  6.00a  contained  5  trajectories. 
The  initial  conditions  were  selected  as  follows:  for  each  value  of  rcut,  a  is  set  to  the  minimum 
value  in  the  range  of  a’s  to  be  investigated  and  8  is  assigned  an  arbitrary  value.  The  particles 
were  placed  in  lattice  sites  of  either  a  bcc  or  fee  crystal  to  generate  the  initial  coordinates.  Initial 
velocities  were  set  to  zero,  and  the  coordinates  slightly  displaced  from  the  equilibrium  lattice 
sites.  The  equations  of  motion  were  integrated  for  50,000  time  steps  at  a  temperature  of  135  K, 
and  a  pressure  of  40  bar.  This  was  found  to  be  sufficient  for  equilibration.  Once  the  system 
equilibrated,  the  equations  of  motion  were  then  integrated  for  an  additional  75,000  time  steps, 
but  s  was  allowed  to  vary  according  to  equation  2.  The  trajectory  was  stopped,  and  the  final 
value  of  s  (denoted  so)  and  the  coordinates  and  velocities  were  used  as  initial  conditions  for  each 
trajectory  with  that  value  of  rcut  over  all  the  values  of  a  and  temperature  studied.  It  is  important 
that  a  suitable  so  is  obtained  so  that  the  system  will  not  prematurely  “vaporize”  (the  volume  of 
the  simulation  box  will  expand  without  bound  during  an  NPT  simulation).  The  problem  is 
circumvented  by  performing  the  equilibration  of  the  system  near  the  boiling  point  (135  K)  and 
using  the  smallest  value  of  a  (in  this  case,  3.35  A)  in  the  range.  This  results  in  the  shortest 
attractive  tail  for  the  interaction  potential  (since  rcut  is  chosen  to  be  a  function  of  a)  over  the 
range  of  a’s. 

Some  care  must  be  taken  not  to  put  the  system  too  close  to  the  boiling  point,  otherwise,  the 
system  might  begin  to  make  the  phase  transition,  which  would  result  in  a  poor  value  for  so.  Each 
trajectory  to  be  used  for  <  s  >  was  integrated  for  25,000  time  steps  beginning  with  the  initial 
conditions  previously  described  and  allowing  s  to  vary  according  to  equation  2.  The  25,000  time 
steps  were  sufficient  for  the  system  to  adjust  to  each  new  value  of  a.  By  this  time,  the  density  of 
the  system  was  oscillating  about  the  desired  value.  After  this  period,  the  instantaneous  density 
was  monitored.  If  the  instantaneous  density  was  within  5  x  10-6  of  the  experimental  value  at  that 
temperature  and  pressure,  the  trajectory  was  stopped,  and  the  value  of  8  at  that  point  was  selected 
as  the  correct  value  for  that  a  value.  Once  this  entire  procedure  has  been  performed  for 
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the  series  of  o’ s  at  a  single  temperature,  the  process  is  repeated  for  that  temperature,  except  the 
initial  conditions  for  the  subsequent  run  (coordinates  and  velocities)  are  those  corresponding  to 
the  immediately  preceding  trajectory.  This  is  repeated  until  the  desired  number  of  values  of  e’s 
that  will  be  used  in  averaging  are  calculated  for  a  given  temperature.  The  a’s  that  were  sampled 
ranged  from  3.35-3.45  A  (in  increments  of  0.025  A)  and  the  temperatures  sampled  ranged  from 
85-145  K  (in  intervals  of  5  K).  Once  the  set  of  e’s  and  a’s  have  been  determined  for  all  of  the 
temperatures,  the  common  set  can  be  obtained  by  superimposing  plots  of  the  curves  showing 
<  s  >  vs.  a  for  all  temperatures  sampled  in  the  study. 

In  order  to  evaluate  the  performance  of  the  modified  Lennard-Jones  potential  using  the 
common  set  obtained  for  rcut  =  6.00a,  time  averages  of  bulk  properties  of  the  liquid  over  the 
temperature  range  85-145  K  were  calculated.  For  each  temperature,  a  trajectory  consisting  of 
250,000  time  steps  was  integrated  to  generate  the  ensemble  averages.  The  initial  coordinates  and 
velocities  of  a  system  of  2,000  particles  were  generated  by  first  setting  up  the  particles  in  an  bcc 
lattice,  and  performing  a  trajectory  integration  for  50,000  time  steps  at  the  desired  temperature 
and  pressure,  using  equation  1  and  rcut  =  6.00a.  The  ensemble  averages  were  obtained  by 
averaging  over  the  calculated  values  for  the  remaining  200,000  time  steps. 

2.2  Monte  Carlo  (MC).  This  technique  can  be  implemented  using  MC  simulations  as  well. 
MC  does  not  require  continuous  forces  at  the  cutoff  distance;  therefore,  this  method  will  be 
demonstrated  using  the  unmodified  Lennard-Jones  potential.  Long-range  corrections  are 
included  to  account  for  all  interactions  of  pairs  whose  distances  exceed  rcut  [30,  31].  The  density 
constraint  can  be  implemented  in  the  same  way  as  was  done  for  MD  using  equation  2.  However, 
the  time  dependence  in  equation  2  is  not  applicable  to  MC,  and  therefore,  it  is  replaced  with  a 
“step”  dependence, 


step 

-step  —  ^0  P  » 

i=l 


(3) 
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where  X  is  the  coupling  strength  as  previously  described,  sstep  is  the  value  for  s  for  that  step  of  the 
Markov  sequence,  so  is  the  initial  value  for  e  and  Ap;  is  the  difference  between  the  desired  and 
instantaneous  bulk  densities. 

Parameter  sets  for  various  values  of  a  over  the  temperature  range  85-145  K  were  determined 
using  NPT-MC  simulations  and  equation  3.  As  in  the  MD  simulations,  these  sets  were 
determined  for  a  values  ranging  from  3.35-3.45  A  in  increments  of  0.025  A'  The  simulation 
box  contained  150  atoms,  placed  randomly  within  the  box,  and  the  size  of  the  box  was  initially 
set  such  that  the  initial  density  was  equal  to  the  experimental  value  at  135  K,  40  bar.  As  in  the 
MD  simulations,  a  is  set  to  the  minimum  value  in  the  range  to  be  studied,  and  s  is  given  an 
arbitrary  value.  The  NPT-MC  simulation  proceeded  for  100,000  moves,  with  s  remaining  fixed. 
This  was  sufficient  to  equilibrate  the  system  to  the  desired  thermodynamic  state.  Imposing  the 
same  temperature  and  pressure  constraints,  the  NPT-MC  simulation  proceeded  for  an  additional 
100,000  Markov  moves  but  the  s  was  now  allowed  to  vary  according  to  equation  3.  The 
sequence  continued  with  a  trial  move  of  the  particle  positions,  which  generated  an  instantaneous 
density.  A  new  value  of  s  was  generated  according  to  equation  3.  The  energy  of  the  system 
using  the  new  value  of  s  was  calculated,  and  the  changes  accepted  or  rejected  according  to  the 
probability  min  [1,  exp  (-W/kT)]  [30, 31]  where 


W  =  P(v_-V0J+(U_-Uj+McTln^,  (4) 

and  P,  V,  and  U  denote  the  pressure,  volume,  and  potential  energy  of  the  system,  respectively. 
The  magnitude  of  the  displacements  of  the  particles  during  the  Markov  sequence  were  chosen 
such  that  50%  of  all  attempted  moves  were  accepted.  If  a  trial  move  was  rejected,  the  properties 
associated  with  the  immediately  preceding  configuration  were  included  for  averaging. 

The  final  value  of  s  at  the  end  of  the  100,000  moves  became  so  and  the  resulting  coordinates 
were  used  for  the  initial  conditions  for  NPT-MC  calculations  to  determine  s  as  a  function  of  a 
over  the  entire  liquid  phase  at  40  bar.  As  in  the  MD  simulations,  this  process  is  repeated  for  the 
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series  of  cr’s  at  each  temperature,  except  that  the  initial  conditions  for  the  subsequent  run 
(coordinates)  are  those  corresponding  to  the  final  configuration  of  the  immediately  preceding 
Markov  sequence.  This  is  repeated  five  times  in  this  study.  The  common  set  is  determined,  as 
before,  by  superimposing  the  plots  of  <  s  >  vs.  c  for  each  temperature. 

In  order  to  assess  the  performance  of  the  Lennard-Jones  potential  using  this  set  of  parameters 
and  those  derived  from  other  studies  [36,  37],  NPT-MC  simulations  of  liquid  argon  at  40  bar 
over  the  temperature  range  85-145  K  were  performed.  The  initial  system  consisted  of  108  atoms 
in  a  cubic  box  where  each  edge  length  was  17.46  A.  Periodic  boundary  conditions  were 
imposed,  and  a  cutoff  distance  equal  to  one-half  of  the  edge  of  the  box  was  used.  Long-range 
corrections  to  the  calculations  were  included.  An  initial  equilibration  simulation  using  the 
Lennard-Jones  parameters  of  Wijker  et  al.  [37]  for  T  =  97  K,  P  =  17  atm  was  performed  in  order 
to  reproduce  one  of  the  values  generated  by  McDonald  and  Singer  [36].  An  excess  of  500,000 
Markov  moves  were  used  for  equilibration,  and  100,000  moves  were  used  in  averaging  the 
results.  The  results  using  this  parameter  set  were  reproduced.  This  coordinate  set  was  then  used 
as  the  initial  set  for  the  calculations  of  properties  at  85  K  using  parameters  obtained  in  this  work, 
from  McDonald  and  Singer  [36],  and  from  Wijker  et  al.  [37],  An  equilibration  Markov  sequence 
of  10,000  steps  was  performed  at  T  =  85  K,  P  =  40  bar  and  thermodynamic  properties  were 
calculated  and  averaged  over  the  next  100,000  steps.  A  new  temperature  was  then  selected,  and 
the  process  repeated  until  the  entire  temperature  range  had  been  sampled.  Longer  Markov  chains 
were  generated,  and  different  sequences  were  initiated,  but  the  results  did  not  vary. 


3.  Results 

Although  argon  has  been  well  described  using  more  complex  interaction  potentials  [38,  39] 
the  purposes  of  this  study  are  (1)  to  illustrate  the  newly  described  method  for  obtaining 
parameter  sets  using  either  MD  or  MC,  and  (2)  to  show  that  the  Lennard-Jones  potential  can 
obtain  accurate  results  over  the  entire  liquid  phase  upon  proper  parameterization. 
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3.1  Molecular  Dynamics  (MD).  Parameters  for  the  modified  Lennard-Jones  potential 
(equation  1)  were  determined  using  the  method  described  in  the  preceding  section  for 
temperatures  of  85-145  K.  The  value  of  <e>/k  and  the  corresponding  a  for  liquid  argon  at  a 
pressure  of  40  bar  at  several  temperatures  is  shown  in  Figure  1.  This  is  presented  for  cutoffs 
rcut’s  of  2.50a  and  6.00a.  Figure  1(a)  shows  <s>/k  vs.  a  using  rcut  =  2.50a,  and  Figure  1(b) 
shows  the  values  using  rcut  =  6.00a.  In  the  figure,  the  common  set  corresponds  where  all  of  the 
curves  intersect  (within  statistical  error).  This  changes  considerably  from  2.50a  until  about 
4.00a,  where  it  slowly  converges  to  its  value  at  6.00a.  The  common  set  at  6.00a  was  used  to 
calculate  all  properties  of  interest  for  the  MD  simulations.  The  error  bars  for  all  points  are  on  the 
order  of  the  size  of  the  symbols  except  for  those  points  at  145  K,  which  show  the  error  bars 
explicitly.  The  deviations  at  145  K  are  due  to  the  proximity  of  the  system  to  the  boiling  point; 
the  system  is  beginning  to  make  the  expected  phase  transition.  Since  the  values  generated  at 
145  K  have  such  greater  uncertainty  than  the  remaining  points,  these  have  not  been  included  in 
the  determination  of  the  common  set.  For  all  values  of  rcut,  the  a  in  the  common  set  is  3.40  A  but 
the  corresponding  e  values  are  very  different.  For  example,  for  rcut  =  2.50a,  the  value  of  <s>/k  is 
155.876  K,  but  the  value  of  <s>/k  resulting  from  the  simulations  in  which  rcut  =  6.00a  is 
1 19.808  K.  The  changes  in  the  value  of  the  s  of  the  common  sets  with  increasing  values  of  rcut 
are  given  in  Table  1,  and  illustrated  in  Figure  2.  The  values  of  s  for  the  common  set  appear  to  be 
converging  as  rcut  increases  (i.e.,  <e>/k  converges  for  a  =  3.40  A  to  «  120  K  as  rcut  is  extended). 
The  error  bars  have  been  eliminated  since  they  are  comparable  with  the  size  of  the  symbols. 


The  performance  of  the  modified  Lennard-Jones  potential  using  the  common  set  obtained  for 
rcut  =  6.00a  was  tested,  in  order  to  keep  the  system  size  small  enough  that  molecular  simulation 
would  not  be  computationally  prohibitive.  NPT-MD  simulations  were  used  to  calculate 
ensemble  averages  for  density,  internal  energy  and  enthalpy  over  the  temperature  range 
85-145  K  at  40  bar.  Table  2  gives  the  calculated  ensemble  averages  for  densities,  internal 
energies  and  enthalpies  along  with  the  experimental  values  [40],  and  the  comparisons  are 
illustrated  in  Figures  3-5,  each  of  which  are  based  on  MC  calculations  using  the  parameters  of 
Wijker  et  al.  [37],  McDonald  and  Singer  [36],  and  those  obtained  in  this  work  along  with  the 
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Figure  1.  Value  of  <  e  >/k  and  Corresponding  o  for  Liquid  Argon  at  a  Pressure  of  40  bar 
for  Several  Temperatures. 


Table  1.  Value  of  <e>/k  for  the  Common  Set  as  a  Function  of  Cutoff  Distance  Using 
Equation  1 


rCui(cr  )a 


_ Z50 _ 

_ 2/75 _ 

3.00 


_ 325 _ 

_  3.50 _ 

_ 400 _ 

5.00 


6.00 


'The  value  of  O  for  the  common  set  is  3.40  A. 


(<S>/k)/K 


155.876 

144.707 

137.652 

133.112 

129.875 

125.437 

121.516 


119.808 


Deviation 


°I 

26 

15 


71 _ 

84 
0.116 
0.069 


0.088 


experimental  values  [40].  Also  shown  in  Figures  3-5  are  the  values  obtained  from  MD 
simulations  using  the  parameters  found  in  this  work.  The  error  bars  have  been  eliminated  since 
they  are  comparable  to  the  size  of  the  symbols.  The  density  (Figure  3),  the  energy  (Figure  4), 
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Figure  2.  Value  of  <e>/k  for  a  =  3.40  A  Plotted  as  a  Function  of  the  Cutoff  Distance  rcut 
(in  cj)  for  Liquid  Argon. 

and  the  enthalpy  (Figure  5)  have  been  multiplied  by  102  for  convenience.  The  calculated  results 
are  within  1%  of  the  experimental  data  for  most  of  the  results  over  the  temperature  range 
investigated.  (Once  again,  the  deviations  at  145  K  result  from  the  system  being  at  the  boiling 
point.)  These  results  are  comparable  with  results  that  were  obtained  with  a  more  complicated 
interaction  potential  [39]. 

The  ability  of  the  modified  Lennard-Jones  potential  using  the  common  set  obtained  at  rcut 
=  6.00a  to  describe  the  self-diffusion  coefficient  has  also  been  evaluated.  In  their  paper  [41], 
Naghizadeh  and  Rice  present  experimental  data  of  self-diffusion  for  liquid  argon  over  a 
temperature  range  of  90-120  K  for  pressures  of  12.9,  57.5,  103.0,  and  135.0  atm.  The  behavior 
of  the  self-diffusion  coefficients  with  temperature  at  each  pressure  was  satisfactorily  described 
using  an  Arrhenius-like  expression.  Naghizadeh  and  Rice  fitted  their  data  to  this  expression,  and 
provided  an  analytical  description  of  the  temperature  dependence  of  the  self-diffusion  coeffcient 
for  liquid  argon  could  be  estimated  at  a  pressure  of  40  bar  for  each  temperature,  to  allow  for  a 
comparison  with  the  MD  simulations  results. 
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Table  2.  Density,  Internal  Energy,  and  Enthalpy  of  Liquid  Argon 
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Figure  3.  Number  Density  of  Liquid  Argon  (in  Atoms  A  3)  as  a  Function  of  Temperature. 


The  self-diffusion  coefficient  from  MD  is  calculated  from  the  mean-square  displacement  and 
is  an  average  over  ten  trajectories  each  consisting  of  20,000  integration  steps.  Table  3  compares 
the  experimental  and  calculated  self-diffusion  coefficients  as  a  function  of  temperature.  The 
calculated  values  have  a  larger  percentage  deviation  (with  a  maximum  error  of  8%)  from 
experimental  values  than  those  of  the  thermodynamic  properties,  but  are  still  in  reasonably  good 
agreement.  However,  Naghizadeh  and  Rice  indicate  that  their  data  contains  uncertainties  of  less 
than  5%,  which  could  account  for  some  of  the  discrepancies. 

3.2  Monte  Carlo  (MC).  Lennard-Jones  parameters  for  liquid  argon  were  obtained  using  the 
method  described  in  section  2  and  NPT-MC  calculations  over  the  temperature  range 
85-145  K,  at  40  bar.  Figure  6  shows  <  s  >/k  over  the  a  range  3.35-3.45  A  in  increments  of 
0.025  A  obtained  from  NPT-MC  calculations  using  this  method.  In  this  figure,  the  common  set 
corresponds  to  a  =  3.40  A  as  it  did  in  the  application  of  the  method  using  NPT-MD.  The  value 
of  e,  averaged  over  all  temperatures  is  <e>/k  =  116.359  ±  0.128  K.  The  value  for  s  is  3% 
smaller  than  that  obtained  by  Wijker  et  al.  [37]  and  is  within  the  error  of  the  value  given  by 
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Figure  4.  Internal  Energy  (in  eV  Atom-1)  as  a  Function  of  Temperature. 


McDonald  and  Singer  (117.2  K  ±  1.4  K)  [36].  The  value  for  a  that  was  calculated  is  close  to 
that  of  Wijker  et  al.  [37]  and  McDonald  and  Singer  (3.40  A)  [36];  however,  when  generating  the 
curves  of  e  as  a  function  of  o  at  each  temperature,  those  for  a  fine  grid  of  o  values  were  not. 
Rather,  these  points  were  generated  at  intervals  of  0.025  A. 


The  parameters  of  Wijker  et  al.  [37]  were  obtained  using  experimental  gas-phase  information 
for  argon,  rather  than  liquid  state  information.  Small  discrepancies  between  experiment  and 
results  from  constant  volume  and  temperature  (NVT)-MC  simulations  using  these  parameters  to 
describe  thermodynamic  properties  in  the  liquid  state  led  McDonald  and  Singer  to  propose  a  new 
set  of  parameters  that  would  better  reproduce  values  in  the  liquid  state  [36].  In  order  to  compare 
the  results  described  herein  with  those  of  experiment  (Wijker  et  al.  [37]  and  McDonald  and 
Singer  [36]),  NPT-MC  simulations  were  performed  to  predict  densities,  internal  energies  and 
enthalpies  for  liquid  argon  at  40  bar,  over  the  temperature  range  85-145  K  using  Lennard- Jones 
potentials  with  the  parameters  of  McDonald  and  Singer  [36],  Wijker  et  al.  [37],  and  the  set  that 
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Figure  5.  Enthalpy  (in  eV  Atom  ')  as  a  Function  of  Temperature. 


Table  3.  Self-Diffusion  Coefficients  (in  cm2-s  ’)  for  Liquid  Argon 


T/K 

DexpH  105 

DmdH  105 

85 

1.681 

1.783 

90 

2.122 

2.281 

95 

2.614 

2.709 

100 

2.954 

3.191 

105 

3.737 

3.764 

110 

4.361 

4.533 

115 

5.022 

5.374 

120 

5.715 

5.950 

15 


130 


Figure  6.  The  Value  of  <  e  >/k  and  Corresponding  a  for  Liquid  Argon  at  a  Pressure  of 
40  bar  for  Several  Temperatures. 


was  derived  using  the  method  described  in  this  work.  The  results  of  all  calculations  are  given  in 
Table  2  for  comparison  with  experiment  and  among  the  predicted  values;  Figures  3-5  provide 
visual  comparisons.  With  the  exception  of  values  at  145  K,  the  calculations  using  the  Wijker  et 
al.  [37]  set  are  in  poorer  agreement  with  experiment  than  potentials  using  the  set  derived  herein 
or  the  McDonald  and  Singer  set  [36],  Properties  calculated  at  145  K  using  both  the  McDonald 
and  Singer  set  and  the  parameters  determined  in  this  work  have  a  significant  deviation  from 
experimental  values.  The  values,  however,  do  not  represent  an  average  over  a  single  phase  since 
the  system  is  actually  fluctuating  between  the  liquid  and  vapor  states  [30],  which  was  evident 
upon  inspection  of  the  instantaneous  densities.  At  temperatures  away  from  this  phase  transition, 
the  properties  predicted  using  either  McDonald  and  Singer  or  our  set  are  comparable,  thus 
affirming  McDonald  and  Singer’s  conclusion  that  the  parameters  generated  by  Wijker  et  al. 
produce  an  attractive  well  too  deep  for  liquid  argon  [36], 
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The  parameters  generated  in  this  work  were  determined  using  experimental  information  of 
the  density  of  liquid  argon  at  40  bar;  in  order  to  assess  the  performance  of  this  potential  at  other 
pressures,  NPT-MC  calculations  were  performed  at  pressures  ranging  from  100-600  atm,  at  100, 
120,  and  140  K  for  comparison  with  experiment.  The  resulting  molar  volumes  are  given  in 
Table  4  and  compared  with  the  experimental  results  of  Streett  and  Staveley  [42].  The  predicted 
values  over  the  entire  pressure  and  temperature  range  are  in  agreement  with  experiment  to 
0.83%. 


Table  4.  Molar  Volume  of  Liquid  Argon  (in  cm3  mol  *) 


|  T/K  | 

P/atm 

100.0 

121 

D.O 

140.0  | 

Experiment 

This  work 

Experiment 

This  work 

Experiment 

This  work 

100.0 

29.69 

29.64 

32.77 

32.88 

37.49 

37.80 

200.0 

29.01 

28.94 

31.52 

31.66 

34.91 

35.11 

300.0 

28.46 

28.35 

30.62 

30.70 

33.20 

33.44 

400.0 

27.98 

27.88 

29.91 

29.93 

32.14 

32.27 

500.0 

27.57 

27.47 

29.33 

29.30 

31.27 

31.35 

600.0 

27.23 

27.08 

28.82 

28.83 

30.58 

30.66 

4.  Conclusion 

It  has  been  shown  that  it  is  possible  to  optimize  a  pair-additive  interaction  potential  using 
molecular  simulations  in  which  the  system  is  constrained  to  reproduce  the  experimental  bulk 
property  over  a  range  of  temperatures  at  a  constant  pressure.  Specifically,  this  was  shown  for 
liquid  argon  using  either  a  modified  or  unmodified  Lennard- Jones  potential,  and  using  the  bulk 
density  of  the  liquid  as  the  constraining  factor.  The  modification  to  the  Lennard-Jones  potential 
used  in  the  NPT-MD  simulations  was  made  to  conserve  the  energy  during  the  trajectory 
integration  and  to  ensure  the  force  was  continuous  at  the  cutoff  distance.  Since  NPT-MC 
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simulations  do  not  require  this,  application  of  the  method  using  NPT-MC  was  performed  using  a 
regular  Lennard-Jones  interaction  potential  with  cutoff. 

NPT-MD  simulations  using  this  method  produced  values  of  <e>/k  =  119.808  K  and 
o  =  3.40  A  with  rcut  of  6.00a  were  found  to  be  the  best  set  of  parameters  over  the  temperature 
range  of  85-145  K.  Using  these  parameters,  excellent  agreement  with  experiment  for  density, 
internal  energy,  enthalpy,  and  self-diffusion  coefficient  was  obtained.  This  suggests  that  density 
is  a  good  constraint.  NPT-MC  simulations  using  this  method  produced  values  of  <  s  >/k 
=  1 16.237  K  and  a  =  3.40  A  for  the  unmodified  Lennard-Jones  potential,  which  also  gave  results 
that  were  in  very  good  agreement  with  experiment  over  the  temperature  range. 

Results  from  the  simulations  using  either  NPT-MD  and  NPT-MC  were  compared  against 
available  experimental  information  and  NPT-MC  calculations  using  parameters  for  the 
Lennard-Jones  potential  suggested  by  McDonald  and  Singer  [36]  and  Wijker  et  al.  [37].  The 
results  using  the  parameters  generated  in  this  study  were  in  better  agreement  with  experiment 
than  those  using  the  Wijker  et  al.  set,  and  were  comparable  to  those  generated  using  the 
McDonald  and  Singer  results.  Additionally,  NPT-MC  calculations  were  performed  using  the 
common  set  of  parameters  determined  for  the  Lennard  Jones  potential  for  temperatures  ranging 
100-140  K  and  pressures  of  100-600  atm  for  comparison  with  experimental  results.  The  molar 
volumes  calculated  using  NPT-MC  were  within  0.83%  of  the  experimental  values.  The  results 
presented  herein  seemed  to  indicate  that  s  and  a  are  independent  of  temperature  and  pressure. 

Although  this  method  should  be  applied  to  other  systems  to  determine  if  it  will  show  the 
same  degree  of  success  in  parameterizing  simple  potential  energy  functions,  the  results  herein 
suggest  that  this  might  be  a  very  useful  and  powerful  tool  for  easily  and  effectively  developing 
model  interaction  potentials  to  describe  real  systems.  Such  models  can  then  be  used  to  obtain 
thermodynamic  and  some  transport  properties  from  computer  simulations. 
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